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We propose a model of nonequilibrium quantum transport of particles and energy in a system con- 
nected to mesoscopic Fermi reservoirs (meso-reservoir) . The meso-reservoirs are in turn thermalized 
to prescribed temperatures and chemical potentials by a simple dissipative mechanism described by 
the Lindblad equation. As an example, we study transport in monoatomic and diatomic chains of 
non-interacting spinless fermions. We show numerically the breakdown of the Onsager reciprocity 
relation due to the dissipative terms of the model. 
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Introduction.- Nonequilibrium systems abound in na- 
ture and still their theoretical description poses numer- 
ous challenges to theory. Nonequilibrium steady states 
(NESS) describe the state of a system maintained out of 
equilibrium by external forces such as e.g., gradients of 
temperature T and chemical potential p, and are char- 
acterized by the emergence of steady flows. One might 
expect that a finite system connected to two infinitely ex- 
tended reservoirs imposing external gradients, will reach 
a NESS after a sufficiently long time [HQ. Nonetheless, 
this is not always the case and the conditions to reach a 
NESS are not well understood in general. Although the 
problem of existence and approach to NESS has been 
discussed since the early days of statistical mechanics, 
rigorous results are limited to a few examples , and 
to near equilibrium regimes. The extension to regimes 
beyond linear response faces the difficulty arising from 
infiniteness of the reservoirs [(| 0] • 

In quantum mechanics the construction of NESS re- 
quires to consider open quantum systems, rendering the 
problem extraordinarily more difficult. The common 
setup is to consider the infinite time limit of the den- 
sity matrix of the finite system S coupled to two infinite 
reservoirs which are in thermal equilibrium at different 
temperatures and chemical potentials, starting from ini- 
tial separable state pl® Ps® PR- One can then study the 
properties of the density operator of the total (infinite) 
system 0, Q , or the reduced density operator for the (fi- 
nite) system S [8(, obtained by tracing out the reservoirs' 
degrees of freedom. The time evolution of the density 
operator is naturally determined by the Von Neumann 
equation. However, dealing with infinite degrees of free- 
dom is in most cases difficult. A second approach based 
on the master equation of the reduced density operator 
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FIG. 1: (Color online) The schematic view of the two- level 
reservoir open system model. 



is more accessible, albeit at the price of several approxi- 
mations such as e.g., Born-Markov (see e.g. Q). 

In this Letter, we propose a conceptually different 
approach for the study of NESS that comprises meso- 
reservoirs with a finite number of degrees of freedom, 
which in turn are thermalized by Markovian macroscopic 
reservoirs using a simple Lindblad equation. This setup 
yields a computationally efficient model in which the 
NESS of the system and of the meso-reservoirs are math- 
ematically accessible, thus allowing to study characteris- 
tic features depending on the number of degrees of free- 
dom of the meso-reservoirs, such as relaxation times of 
the system and of the meso-reservoirs, or correlations 
among them. It is interesting remarking that finite meso- 
reservoirs in principle allow the creation of entanglement 
among the different degrees of freedom of the system 
solely through its coupling with the meso-reservoir [loj ]. 
thus making our setup much richer. 

The model.- We consider a one-dimensional quantum 
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chain of spinless fermions coupled at its boundaries to 
meso-reservoirs comprising a finite number of spinless 
fermions with wave number k (k G {1, . . . , K}). The 
Hamiltonian of the system can be written as H = Hg + 
H L +H R + V, where 
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and where {£./} are the nearest neighbor hoppings, U is 
the onsite potential, v is the coupling between the sys- 
tem and the meso-reservoirs, and Cj,cj is the annihila- 
tion/creation operator for the spinless fermions of the 
chain, while ak.a, &\ a are that of the left and right meso- 
reservoirs. 

The key idea of our model is to enforce the finite meso- 
reservoirs to equilibrium using simple Lindblad dissipa- 
tors [Ill4l3j ). If the couplings to the Lindblad dissipa- 
tors are small, then we can interpret these terms as com- 
ing from tracing out infinitely extended (super-)reservoirs 
(schematically depicted in Fig.[TJ. Thus, the density ma- 
trix of the total setup evolves according to the many- 
body Lindblad equation 



dt 



P = -i[H,p] 



(2) 



^ \^ L k,a,mpLl <a!m - {Ll a m L k , a ,m, p} 
k,a,m 



Fk,a.l 



\A\M a l!a, Lk, a ,2 — \J^k 
7(1 - F a (e k )), r k , a: 2 = jF^Ck) , 



2 a fca 



where Lk tOL ,i and Lk, a ,2 are operators representing the 
coupling to the super-reservoirs, F a (e) = (e^ 6- ' 1 ") + 
are Fermi distributions, with inverse tempera- 
tures f5 a and chemical potentials [i a , and [•,•] and {•,•} 
denote the commutator and anti-commutator, respec- 
tively. Generalization to spinfull fermions (e.g. elec- 
trons) is straightforward. The parameter 7 determines 
the strength of the coupling to the super-reservoirs and, 
as we shall show later, needs to be fine-tuned in order 
to ensure the applicability of the model. We note that 
the parameter 7 determining the rate of relaxation of the 
meso-reservoirs to equilibrium, in general can depend on 
temperature, chemical potential, and wave-number k of 
the meso-reservoir modes. Here we consider the simplest 
model possible and take 7 as constant. We stress that our 
model does not rely on the usual weak-coupling assump- 
tion needed for the physical derivation of the Lindblad 
master equation Q, thus 7 does not need to be a small 
parameter. 



Results.- We have studied monoatomic (tj = t) and 
diatomic (t2j-i = ^A, t 2 j = ^b) chains. Unless spec- 
ified differently, we have set 7 = 0.1, v — 0.03, ei = 
-20, e K = 20, K = 200, t = 3, and [i L = -pL R = P- 

When the system and the meso-reservoirs are decou- 
pled (v = 0), each non-interacting mode of the meso- 
reservoir is thermalized separately with the prescribed 
Fermi-Dirac occupation number 14|. For small v, »<<], 
we thus expect the distribution of occupations in the 
meso-reservoirs to be close to Fermi-Dirac, as shown in 
Fig. [21a). A remarkable feature of our model is that we 
can monitor and control the difference of occupation dis- 
tributions {a\ a a,ka) — F a (ek) by changing the coupling 
parameters » or 7. Here and bellow (•} represents an 
expectation value with respect to the NESS. 

Looking at the occupation density inside the system we 
observe that (c|cj) is almost a constant, i.e. independent 
of j, except at the edge of the chain, indicating that 
ballistic transport is expected. 

We now discuss the relaxation times of the system. As 
shown in [HI], the spectrum of the evolution superop- 
erator is given in terms of the eigenvalues f3j (so-called 
rapidities) of the matrix X: 

Okxu Oa"xA"\ 
0», v „ nX K I (g) E2 
KxK v Kxn 

where Oi X j and Ej denote i x j zero matrix and j x j 
unit matrix, a y is the Pauli matrix, and H is a ma- 
trix which defines the quadratic form of the Hamilto- 
nian, as H = d T Hd in terms of fermionic operators 
d T = {a 1L , ■ ••o Jti ,Ci, • • • ,c n ,a 1R , ■ ■ -a KR }. Fig. HJb) 
shows a typical spectrum for the monoatomic chain. In- 
terestingly, there is a clear separation of the relaxation 
times into slow and fast normal modes. The number of 
slow modes n s , with eigenvectors localized in the system 
part, is n s « 2n, and the number of fast modes nj « 
with eigenvectors localized in the meso-reservoirs. 

We have not observed rapidities with zero real part, 
and thus according to Ref. 1J, [l5| , there exists a unique 
steady state for the range of parameters used in this Let- 
ter. Next, we study the NESS averages (■} of obscrvables. 
The quadratic obscrvables are given in terms of the so- 
lution of the Lyapunov equation [il^ : 

(wjWk) = Sj, k - 4iZ jtk 

w 2j -i = c 3 + ct, w 2j = i(cj - ct) 

X T Z + ZX = M, 




M, = --diag{r iL , • ■ ■ ,T KL , lx „, T 1R , ■ 



, r KR } g <Ty 



whereas Wick's theorem can be used to obtain expecta- 
tions of higher-order observables. Let us now discuss the 
deviation of occupation numbers from the Fermi distribu- 
tion in the meso-reservoirs. Lindblad equation of motion 
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FIG. 2: (Color online) Panel (a) shows (mode) fc-dependence 
of the deviation from the Fermi function (a ka a ka ) — F a (tk), 
(red/blue curve is for the left/right meso-reservoir, and verti- 
cal lines indicate the energy band, see main text for parame- 
ter values), and Panel (b) shows the rapidity spectrum for the 
monoatomic chain indicating the separation of time scales. 



results in the following identity, 2"f |(a£afc) — F a (ek) \ = 

(J%), where J£ is the fc-th level contribution to the cur- 
rent from the a-meso-reservoir to the system, J£ = 



iv(aU 



ci, c R 



This relation 



can be used to determine the Lindblad dissipator from 
physical obscrvablcs such as currents and occupation 
numbers. It is interesting to remark that this rela- 
tion can be interpreted as a type of Landauer formula, 
i-e-, (Jfot) = J2k 2 ^{( a kL a kL) - F L (e k )}, where jf ot is 
the particle current from the left reservoir to the sys- 
tem, and shows explicitly how the nonequilibrium sit- 
uation modifies the Fermi distributions. A similar ex- 
pression holds for the current from the right reservoir 
to the system. It follows that although the distribu- 
tion functions (a£ cifc Q ) are slightly modified by the pres- 
ence of the system and the other reservoir, the inte- 
grated difference from the Fermi distribution summed 
over both meso-reservoirs satisfies the charge conserva- 
tion, i.e., J2k, a {( a L a fca) - F a (e k )} =0@. 

We now turn our attention to the particle and energy 
currents in the system. The particle current is defined 
through the conservation law of number of particles: 



dc] Cj 



J? - J?, (2 < j < n - 1) 



dt 3 
Jj = itj{c]c j+1 -c] +1 Cj), (l<j<n-l) 

and the energy current is defined through the conserva- 
tion law of local energy: 

■>';■ (2<i<»-i) 

J^^ilH^Hj], (2<j<n) 

Hj = tjC\Cj + i + tjc' j+1 Cj +Uc]cj, (1 < j < n) 
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FIG. 3: (Color online) Panel (a) shows the K dependence 
of the particle current for the monoatomic chains (red: 7 
f 1. ()(;() : , green: 7 = 0.001, blue: 7 = 0.01, pink: 7 = 0.1, cyan: 
7=1, yellow: 7 = 10), and panel (b) shows the rapidities 
for 7 = 0.01 (Using the symmetry with respect to real axis, 
we plot rapidities in positive/negative imaginary plane for 
K = 216/K = 220). Those two parameters are indicated 
in figure (a) as red (K = 216) and blue (K — 220) circles. 
Panel (c) shows the 7 dependence of the particle current for 
the monoatomic chains. Panel (d) shows fi dependence of the 
particle current. The red line is for the monoatomic chains 
(t — 3), and the blue line is for the diatomic chains (tA = 
3, t B = 6). 



where, by definition, c n+1 = 0. We have checked 
that both the particle and the energy current converge 
to a constant by increasing the system size n for the 
monoatomic chains. For the diatomic chains, currents 
converge to different constants, depending on the par- 
ity of integer n. Therefore, ballistic transport is indeed 
achieved. 

Fig. [3]Ja) shows the K dependence of the particle cur- 
rent for the monoatomic chains with 7 ranging from 
0.0001 to 10. The particle current increases linearly in 
the number of meso-reservoir modes K, with notable fluc- 
tuations in K observed for some range of 7 (which will 
be discussed later). We see that the particle current is 
roughly 7-indcpendent for 7 6 [0.001, 1], and is somewhat 
smaller for very small or very large 7. The dependence 
on 7 of the particle current, shown in Fig. [3Jc), exhibits 
a plateau starting at 7 w v, and one should take 7 in 
the plateau region [v, \\H\\] to have a physically mean- 
ingful model. Fig. [3]Jd) shows the dependence of the 
particle current on the bias of the chemical potential (/i) . 
We find, as expected, initially linear increase in fj, which 
slows down for larger fi. 

The fluctuations in the current as a function of K oc- 
cur because for some values of K the coupling between 



4 




For sufficiently small thermal and chemical gradients, 
the particle and heat current defined as Jq = Jg — p,Jp 
(p = + u r)/ 2), depend linearly on the external gra- 



dients as 17Hl9j 



FIG. 4: (Color online) Panel (a) shows the tp dependence 
of the particle current for the chains. The red line is for 
Pl = /3r = 1, and fi — 1, and the blue line is for /3l = 
0.5, f3p = 4, and /i = 0. Panel (b) shows tp dependence of 
the Onsager coefficients (color code: Lqp - red, Lpp - black, 
Lqq - green, Lpq - blue), and (c) of ZT (comparing K — 200 
case (green) to K — 100 (red) and K = 500 (blue)), for the 
diatomic chains (tA = 3). Chains have exponentially small 
current outside it £ [1,7] indicated in the figure. Panel (d) 
shows the 7 dependence of asymmetry Lqp /Lpq — 1 for the 
diatomic chain with tA = 3, tp — 6 (dashed line indicates 
linear growth). Note that asymmetry grows linearly with 7 
in the physical regime 7 > v. 



the meso-reservoirs and system is stronger. Intuitively 
this happens when an eigenenergy of the meso-reservoir 
coincides with one of the system, resulting with some 
eigenvector of the matrix X being supported both in the 
system and the meso-reservoirs, breaking the separation 
into slow and fast modes. When the density of states 
of the meso-reservoir becomes dense enough such that it 
is smooth (average level spacing smaller than the width 
given by the coupling 7 to the super-reservoirs) , the cur- 
rent becomes a smooth function of K or equivalently of 
the density of states. These ideas are confirmed by the 
numerical results in Fig. [3ta,b). We note however, the 
current fluctuations are strong only in the 'non-physical' 
regime 7 < v. 

Fig. E|a) shows the dependence of the particle cur- 
rent driven by either thermal or chemical gradients, on 
the hopping strength (tp) for diatomic chains. Where 
the particle current is negative, the energy band of the 
chain is located entirely below the fermi energy of the 
cooler reservoir, yielding a particle current flow from the 
cold to the hot reservoir. Such crossed transport (and its 
counterpart of the energy current driven by the chemi- 
cal gradient), can be exploited to pump heat or particles 

El. 



J Q = i QQ A/3 + iQpA(-/3/i) , 
J P = LpqA/3 + L PP A(-^) , 

where A/3 = f} R - f3 L and A(fin) = PrUr - Pl^l- 
The second law of thermodynamics imposes definitc- 
positiveness of the matrix of Onsager coefficients L, 
which implies Lqq > and Lpp > 0, and if the dy- 
namics is time-reversible, the Onsager's reciprocity rela- 
tion Lpq = Lqp holds. In Fig. [4] we consider diatomic 
chains with tA = 3 and show the dependence of vari- 
ous properties of L on the other hopping parameter tp- 
Fig. Hljb) shows the tp dependence of all Onsager coef- 
ficients, whereas in Fig. |4tc) shows the thermoelectric 
figure-of-merit ZT = LpQLQp/detL (20[. One sees that 
while non-diagonal elements Lqp or Lpq can be some- 
times negative, Lqq, Lpp are always positive. Moreover, 
ZT reaches large values only for disproportionate hop- 
ping rates tA and tp- 

On the other hand, systems interacting with envi- 
ronments inevitably include irreversible processes, which 
breaking down time-reversibility invariancc and thus, the 



validity of the Onsager reciprocity relation 21 1. Fig.|4|d) 
shows the 7 dependence of Lpq / Lqp, and we see that the 
relation is roughly linearly broken by increasing 7, and 
thus, we conclude that the Onsager reciprocity relation is 
satisfied only if there is a time reversible dynamics for the 
total system, including the super-reservoirs, which is the 
case in which the Onsager reciprocity relation is rigor- 
ously proved 22( . We recall that 7 <C 1 is one of the nec- 
essary conditions for deriving the Lindblad equation, by 
taking a partial trace of the unitary time evolution, thus 
\LpqI 'Lqp — 1| can be understood as an error- indicator 
due to the weak-coupling assumption. We remark that 
the transport coefficients for the diatomic chains (shown 
in Fig. [3jd) and Fig. [4]) are non-smooth functions of tp 
at the mono-atomic point tp = tA- 

Conclusion.- We have introduced a model in the 
framework of open quantum systems with finite meso- 
reservoirs which nicely describes non-equilibrium steady 
states of quantum chains. We have checked that this 
model has a regime which is robust with respect to the 
strength of the Lindblad dissipators, and the occupation 
number distributions of meso-reservoirs are close to the 
Fermi distribution, the difference being determined by 
the particle current. We observed that the decay times 
of normal modes show a clear separation into slow and 
fast decaying modes. In a certain regime we find a pos- 
sibility of strongly fluctuating currents (as a function of 
any generic parameter) which is attributed to existence 
of well separated intermediate (resonant) decay modes. 
We find possibility of negative non-diagonal elements of 
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the Onsagcr matrix but confirm the positivity of the full 
Onsagcr matrix. The Onsager reciprocal relation was 
shown to be correct only for weak coupling 7 <C 1, and 
the symmetry is broken linearly as a function of 7. 
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